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X-ray Compton scattering tomography 


By James Webber 


Abstract 

We lay the foundations for a new fast method to reconstruct the electron 
density in x-ray scanning applications using measurements in the dark field. 
This approach is applied to a type of machine configuration with fixed energy 
sensitive (or resolving) detectors, and where the X-ray source is polychromatic. 
We consider the case where the measurements in the dark field are dominated by 
the Compton scattering process. This leads us to a 2D inverse problem where we 
aim to reconstruct an electron density slice from its integrals over discs whose 
boundaries intersect the given source point. We show that a unique solution 
exists for smooth densities compactly supported on an annulus centred at the 
source point. 

Using Sobolev space estimates we determine a measure for the ill posedness 
of our problem based on the criterion given by Natterer in [12]. In addition, 
with a combination of our method and the more common attenuation coefficient 
reconstruction, we show under certain assumptions that the atomic number of 
the target is uniquely determined. 

We test our method on simulated data sets with varying levels of added 


pseudo random noise. 


1 Introduction 


In this paper we investigate the potential for the use of incoherent scattered data for 
2D reconstruction in x-ray scanning applications. The use of scattered data for image 
reconstruction is considered in the literature, typically for applications in gamma ray 
imaging, where the photon source is monochromatic BEIEI. However, in many 
applications (e.g security screening of baggage) a type of x-ray tube is often used 
that generates a polychromatic spectrum of initial photon energies (see section [^for 
an example spectrum). There has been recent interest in the use of energy sensitive 
detectors in tomography m, and in the present paper their application is key to the 
ideas presented. 

Our main goal is to show that the electron density may be reconstructed analyt¬ 
ically using the incoherent scattered data and to lay the foundations for a practical 
reconstruction method based on our theory. We apply our method to a machine con¬ 
figuration commonly used in x-ray CT. In addition, by use of the reconstructed density 
values in conjunction with an attenuation coefficient reconstruction, we show under 
the right assumptions that the atomic number of the target is uniquely determined. 

For a photon incident upon an electron Compton (incoherently) scattering at an 
angle u with initial energy Ex, the scattered energy Eg is given by the equation: 

rp _ __ /-I N 

l + (EA/i?o)(l-cosen) ^ ^ 

where Eq ^ SllkeV is the electron rest energy. Equation ([^ implies that cn remains 
fixed for any given Eg and Ex. So in the case of a monochromatic source, assuming 
only single scatter events, for every fixed measured energy Eg (possible to measure 
if the detectors are energy-resolved) the locus of scattering points is a circular arc 
intersecting the source and detector in question. For example, refer to niEi. 

In an x-ray tube a cathode is negatively charged and electrons are accelerated by a 
large voltage (E^ax kV) towards a positively charged target material (e.g Tungsten). A 
small proportion of the initial electron energy 1%) is converted to produce photons. 
Due to energy conservation, the resulting photon energies are no more than Emax keV. 
So in the polychromatic source case, again assuming only single scatter events, for 
each given data set (photon intensity recorded with energy Eg), the set of scatterers 
lie on a collection of circular arcs intersecting the source and detector points. Together 
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these form a toric section in which the photons scatter, with a maximum scattering 


angle Wmax given by: 


See hgure [T] below: 


cos (cUn 


= 1 


Eq (-Emax — Eg) 
E.E, 


s-^max 



( 2 ) 


Figure 1: A toric section T in which the photons scatter with tips at source and 
detector points s and d. 

In the present paper we consider a setup consisting of a ring of hxed energy sensitive 
detectors and a single rotating fan beam polychromatic source. See hgure With this 
setup we can measure photon intensity in the dark held. We image an electron density 
/ : —)■ M compactly supported within the detector ring (the blue and green circle 
in hgure [^, with / > 0. If we assume an equal scattering probability throughout the 
region R = D H supp (/) leaving only the electron density to vary, and if we assume 
that the majority of scattering events occur within R, then in this case the integral of / 
over D is approximately determined by the scattered intensity recorded at the detector 
d with some hxed energy Eg. See the appendix for an example application where these 
approximations are valid. With these assumptions and with suitable restrictions on 
the support of /, we aim to reconstruct / from its integrals over discs whose boundaries 
intersect a hxed point, namely the source at a given position along its scanning path. 

In section we present a disc transform and go on to prove our main theorem 
(Theorem [^, which explains the relationship between our transform and the straight 
line Radon transform. As a corollary to this theorem, with known results on the Radon 
transform, we show that a unique solution exists on the domain of smooth functions 
compactly supported on an annulus centred at the origin. Here based on the criterion 
of Natterer in na and using the theory of Sobolev space estimates, we determine a 
measure for the ill posedness of our problem. 
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light field 



Figure 2: An example machine configuration is displayed. A disc D whose boundary 
intersects a point source s and a detector d forms the scattering region R. The source 
s travels along the circular path shown. Detectors under direct exposure to the initial 
x-ray beam are said to be in the light field. Detectors not in the light field are said to 
be in the dark field. 

In section we discuss a possible means to approximate the physical processes 
such as to allow for the proposed reconstruction method. Here we also present a least 
squares fit for the total cross section (scattering plus absorbtion) in terms of Z (the 
atomic number). From this, we show that Z is uniquely determined by the attenuation 
coefficient and electron density. 

In section we apply our reconstruction formulae to simulated data sets, with 
varying levels of added pseudo random noise. This is applied to the given machine 
configuration. We recover a simple water bottle cross section image (a circular region 
of uniform density 1) and reconstruct the atomic number in each case using the curve 
fit presented in section To give an example reconstruction of a target not of uniform 
density, we also present reconstructions of a simulated hollow tube cross section. 
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2 A disc transform 


In this section we aim to recover a smooth function compactly supported on an annulus 
centred at the origin O from its integrals over discs whose boundaries intersect O (the 
given source position). 

Let denote the set of points on the disc whose boundary intersects the origin, 
with centre given in polar coordinates as (p/2,0). See figure]^ Let (hi) be the 
set of smooth functions on hi C M" and let (f2) denote the set of smooth functions 
compactly supported on Let Z~^ = x and for a function in the plane 
/ : —)■ M, let F : —)■ M be defined as F{p,6) = / (pcos 0, psin 0). Then we 

define the disc transform Vi : (M^) —)■ C°° (Z~^) as: 

TT COS 9 

'^if(p,(/>)= II fdA= I I pF{p,9 + (j))dpd9 (3) 

JjDr J-iJo 

p 



Figure 3: A disc Dp^p with its boundary intersecting O. 
After making the change of variables: 

p = rcos0, 9 = xIj, dpd9 = cos'ipdrd'ip 

in equation ([^, we have: 


'Dif{p,<P) = 


r cos^ tjjF (r cos 0,0 + 0) d0dr 


(4) 


(5) 


We now present further definitions which will be important in the following subsection 
(section [2T| ), where we provide our Sobolev space estimates. Let Z = M x denote 
the unit cylinder in Then we define T >2 : (M^) —)■ C°° (Z) as follows: 

(P, 0) P > 0 

F2/(p,0) = <! p = o (6) 

(-P, 0 + tt) p<0 
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which is piecewise continuous as a function of p. We can remove this discontinuity by 
adding the function: 

c (0) p > 0 


c (0) sgn (p) 


0 p = 0 


(7) 


—c (0) p < 0 

where c (0) = Pi/(o , 0 ) ^ define V : (7“ (M^) —)■ C°° (Z) as: 


{P, 0) = 7 ^ 2 / (P, 0) + c (0) Sgn (p) 


( 8 ) 


Let Lp,^ = {(x, p) G : X cos 0 + p sin 0 = p} be the set of points on a line. Then we 
dehne the Radon transform R : (M^) —>■ C°° (Z) as: 

Rf (p, 0) = / fds (9) 

We are now in a position to prove our main theorem, where we give the explicit relation 
between R and the Radon transform R for smooth functions on an annulus. 


Theorem 1 . Let Ar^^r2 = G : ri < |x| < r2} be the annulus centred on O with 
inner radius ri > 0 and outer radius r 2 . Let f G (^i,r) for some r > 1 and let 
f eC^ {Ai/r,i) be defined as /(x) = j^f . Then A-Vf = -Rf. 

Proof. Let F and F be dehned as F (p, 0) = f {p cos 0, p sin 9) and F (p, 9) = f {p cos 9, p sin 9). 
Then from our definition of /, we have F (p, 9) = -^F I ^,9). Now we have: 


d 1 


n I cos'll! , , , , , 

cos"^ 'ifF -, 0 + 0 d0 


p 


= -p 


p 


-,0 + 0 


d0 

^ cos 0’ ^ ^ J COS^ 0 

= —Rf (p, 0) for p > 0 


( 10 ) 


and hence §RAf (0, 0) = -§fL)if (0,0) = (0, 0 + vr). So the partial derivative of 


dp 


dp 


Vf with respect to p exists and is continuous for all p G M, and ^Rf = —Rf- 


□ 


We now aim to prove injectivity of the disc transform R on the domain of smooth 
functions compactly supported on an annulus. First we state Helgason’s support 
theorem [6]. 
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Theorem 2. Let X be a compact convex set in M" and let f be continuous on 
If Rf = 0 for all p and cf such that fl X = 0 and f is rapidly decreasing, in the 
sense that: 

\x\^f{x)^0 as |x| —)■ oo Vfc G N (11) 

then f (x) = 0 for all x ^ X. 


Corollary 1. Let / G (^i,r) for some r > 1, and let Zr = {{p, (f) E Z ■. 1/r < p < 
1, (/) G [0, 27r]}. Then / is nniqnely determined by Rf known for all {p, (f) G Z^.. 


Proof. Let: 

/ G {/ G Co^ {A,,r) ■.Vf = 0 for all {p, 0) G (12) 

and let / be defined as in TheoremThen by Theorem]^ we have: 

/ G {/ G {Ayr,i) ■.Rf = 0 for all {p, 0) G (13) 

and hence / is rapidly decreasing. Let X = {x G : |a;| < l/r}. Then X is clearly 


compact and convex. By (13), Rf = 0 for all p and 0 such that fl X = 0. So 
by Helgason’s support theorem, we have that f{x) = 0 for all x 0 X. The result 
follows. □ 


For the proposed machine configuration, we can define the set of points within the 
detector ring formally as Dr = {{x,y) G : x‘^ + {y — {r + 1) /Zf < (r — 1)^/4}, 
where r > 1 depends on the machine specifications (i.e the detector ring radius and 
the source path radius). We now have: 


Corollary 2. Let / G (7“ (Hr)- Let dDr denote the boundary of Dr, and let = 
D\ ,L\ Dr- Then the values of D\f for p and 0 such that: 

p ’V 

Rp,0^0 and ODi^^ndDr (14) 

determine / uniquely. 


Proof. We consider two cases. If Dr C Di ^^ then Vif (p, 0) = Vif 7r/2), which is 
known as condition (14) is satisfied for p = 1/r and 0 = 7r/2. If Dr ft Di ^ = 0, then 
Dif [p, 0) = 0. In any other case, Dif is known by our assumption. Hence we have a 
full data set for Dif and hence for Df. The result follows from Corollary □ 
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So for the proposed application, we see from the above corollaries that for any 
given source position, the incoherent scattered data is sufficient to reconstruct the 
target density uniquely. 

2.1 Sobolev space estimates 

In this section we provide Sobolev space estimates for the disc transform V. From 
these we obtain an upper bound for the least squares error in our solution / in terms 
of e, where e is an upper bound for the least squares error in our measurements. First 
we dehne our Sobolev spaces and the norms which will be used in our estimates. 

Let C M” be an arbitrary domain and let (fl) denote the set of square in- 
tegrable functions on fl. We dehne the Fourier transform of a function / G (M"') 
as: 

/(O = (2vr)“"''^ / /(a;)e"“'«dx (15) 

Jm." 

Then we can dehne Sobolev spaces if" [W^) of real degree a G M as: 

ii" (M*^) = {tempered distributions / : (l + / (0 e (16) 

with the norm: 

ii/irH-d...) = / (i+i?n"i/«)rdf (i?) 

For functions on the cylinder Z C we have the norm: 

\\f\\H^{z)= [ + (18) 

Jr 

where the Fourier transform of Rf is taken with respect to the variable p G M. We 
now state some preliminary results on Sobolev spaces and the Radon transform which 
will be used in our theorems m pages 11 and 203]. 

Theorem 3. For f E S (M^), where S (M^) is the Schwartz space on we have: 

= a eR (19) 

where <F = (cos 0, sin 0) and the Fourier transform of Rf is taken with respect to the 
p variable only. 


Theorem 4. Let k = (/ci,..., kn) he some multi index and let 

dx-^ 

where the ^ are defined in the weak sense. Let m be a positive integer and let a G 
(0,1). Then for a = m + a, the norm 0 is equivalent to the norm: 

ll/ll?,.,= ll/ll?,.,+ E // (20) 

when f2 = M"'. 


We now prove a slice theorem for the disc transform V. 


Lemma 1. Let f G (^i,r) for some r > 1 and let f be defined as in Theorem^ 
Then we have: 

— icrP/(cr, 0) = / ((7$) , a G M (21) 

where <h = (cos 0, sin 0) and the Fourier transform of T>f is taken with respect to the 
p variable. 


Proof. Let 1)2 and c (0) be as dehned in section Then we have: 


Rf (c^, 0) = (27r) 


- 1/2 


i?/(p,0) 


rO poo 

/ Rfip,fi)e-^P^dp+ / i?/(p,0)e-*^"dp 

'-0O Jo 

I ^I>2/ (P, <P) e-’^dp + r (P- 

• poo 

^ - (2r)-‘'^ (Vif (0-,4>) - P2f (0+,<:’)) - I V2f{p,^)e-”dp 


= -(27r)-‘''^ 


(27r 


(27r) 

= [ . + ^ 2 / 0) ] 

za [zn) ' 

= -i(^Rf (o-.0) 


The resnlt follows from the Fourier slice theorem. 


( 22 ) 

□ 


In [T^ page 92] Natterer explains why it is reasonable to consider picture densities 
as functions / of compact support in iL“ (M”) with a < 1/2. He then gives a bound 
for the least squares error in his reconstruction from plane integral data in terms of p, 
where ||/|| hc « < p. With this in mind we will show that the map / —>■ / is bounded 
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and has a bounded inverse from iJ" —)■ for 0 < a < 1. First, from [T^ page 204], 

we have the lemma: 


Lemma 2. Let x ^ (M"') and let f G (ML). Then the map f ^ xf 'Is bounded 

in for any a G M. 

Now we have our result: 

Lemma 3. Let Dr be as defined in section^ Let f G (Dr) for some r > 1 and let 
f be defined as in Theorem Then there exist constants c (a) and C (a) such that: 

c(«) Il/||j7“ < \\f\\H°‘ < C (a) ||/||h“ (23) 


for any 0 < a < 1. 


Proof. Let XOr ^ Cff (M^) be 1 on Dr and let x = XDr\A‘^- Then by Lemmawe 
have: 


Cl (a) 


> \\XDr\x\^f\\l. 


= llklVll 




2 -L 
L2 + 


\x\^f{x) - blV (y) 


-dxdy 


> WfWh + ( 1 ^^) 
>C2(«) WfWl^ 


X — l/|2+2a 

2N2a-2 ff \f(x)-f(y) 


(24) 


k - y\ 


2+2a 


-dxdy 


for 0 < a < 1, This proves the left hand inequality. The right hand inequality can be 
proven in a similar way. □ 


Before proving the main theorem of this section we state the interpolation inequal¬ 
ity for Sobolev spaces on M"' [121 page 203]. 

Lemma 4. Let f G H'^ (M*^). Then we have: 

1I/IIh.(«») < Il/llf^,»)ll/ll5^r, (25) 

for any a < ■j < (3. 


Now we have our main theorem for this section: 
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Theorem 5. Let f G {Dr) for some r > 1. Then we have: 


3/2 


0 




( 26 ) 



for any 0 < (3 < 1 with ||/|| h /3 < P- 

Proof. Let / be defined as in Theorem and let $ = (cos 0, sin 0). Then we have: 

2p/lli.a-i,i]x5i)=2 [ f \Vf{p,cl>)\Mpdcf 

Js^ J-1 

~Is^I 

poo 

\Df {a, 0) + Vf {-a, 0) |Mcrd0 

J —oo 

= -/(-^'^) rdad0 

Js^ J -oo ^ 

>2ti f f (l + ^ 1/((j$) — / (—(j$) |^dcrd0 

JJ —oo 

= 4:71 f f (l + ^ 1/((j$) — / (—(j$) |^d(jd0 

isi Jo 

After making the snbstitntion ^ = a4>, we have: 

^iiwiii.,[-uixs.) > i ler' (1 + 1/(0 - /(-o Pd/ 


P-||2 


> 11 /-/ »„-i 


(27) 


(28) 


where / {x) = f {—x). Applying the interpolation ineqnality with a = —3/2 and 
7 = 0, yields: 


2||/||l2(]r2) < 2||/||/,2(]r2) 


for any 0 < 0 < 1 with 


= 11 /-/ 

/3 _ 3/2 

<ii/-/-ii:xji/-/-iisrM- 


L2(u2) since / G {Dr) 

13 _ _ 3/2 

^3/2 
3 

H-'ii 


) 


3/2 ^ 3/2 

/3 3/2 


0 3/2! 


(29) 


□ 


Corollary 3. Let / G C“ (D^) for some r > 1 and let g = P/. Let G ([—1,1] x 
be snch that \\g‘^ — fi'||L2([_i_i]xsi) < e. Then for any /i,/2 G C“ {Dr) which satisfy 
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W'^f - fi'iL2([_i,i]xsi) < e, we have: 


3/2 p 

ll/l -/2i|L2(K2) < c(/3)p/5+3/2eW (30) 

for any 0 < /3 < 1 with < p for i = 1, 2. 

We can interpret this last corollary to mean that given some erroneons data g‘^ 
which differs in the least sqnares sense from Vf absolntely by e, the least sqnares 

3/2 p 

error in onr solntion is bonnded above by c (/d) for some constant c{f3) 

with the a-priori knowledge that i|/||_ff/3 < p. 

In [12] Natterer nses the valne /?/ (a -f /I) as a measnre for the ill posedness of his 
problem and gives his criteria for a linear inverse problem to be modestly, mildly or 
severely ill posed. If we set (3 close to 1/2, then based on these criteria the above 
argnments wonld snggest that onr problem is mildly ill posed, bnt more ill posed than 
the inverse Radon transform, which we wonld expect given that the disc transform V 
is a degree smoother than R. 



Figure 4: A representation of the detector ring dD^. in the proposed coordinate system. 
The polar angle 9j e [0, 27r] determines the detector position dj. 

Another source of error in our solution can be due to limited sampling of the data. 
In practice the number of detectors will be hnite. Let us parameterize the set of points 
on the detector ring dDr = {(a/, y) G + (p — (r -|- 1) /2)^ = (r — 1)^ /4} in terms 

of a polar angle 0, and let the hnite set of polar angles 0 = {9i,... ,9n} determine 
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a finite set of detector positions {di,... ,dn} G dDr. See figure Then for every 
(j) G [0, 27r] we can sample Vf{p, 0) for: 


P = Pj 


r cos 6j sin 0 + (1 + r sin 6j) sin 0 
+ 1 + 2r sin 6j 


1 < j < n 


(31) 


where pj is such that {H(r — 1) cos{9j), (1 + sin 9j)r + 1 — sin 9j)} C dDj_ .ndD^ for 
1 < j < n. 

From [121 pages 204 and 42] we have: 


Lemma 5. Let 12 C M" 6 e hounded and sujjiciently regular. For h > 0 let 12^ be a 
finite subset of 12 such that d (12,12^) < h, where d is the is the Hausdorff distance 
metric between sets. Let a > n/2 where a = m + a for some integer m and 0 < a < 1, 
and for f G (12) define the seminorm: 


I/I 


2 



D^f {x)-D^f{y) 


x-y\ 


n+2(T 


2 

-dxdy 


(32) 


Then there is a constant c such that: 


WfWi'^in) < (33) 

for every f G (12) which is zero on 12^. 

Theorem 6. Let 12„ he the unit ball in M"'. For every a there exist positive constants 
c (cK, n) and C {a, n) such that for f G (7“ (12„) 

c{a,n) \\f\\H<^(n„) < \\Rf\\Ho.+(r^-i)/2^z) <C(a,n) \\f\\H-in„) (34) 

From these we have the theorem: 

Theorem 7. For each G [0, 27r] let Iff, C [—1,1] be a finite subset of the unit interval. 
Let: 

h = sup d{I^, [-1,1]) (35) 

<t> 

where d is the Hausdorff distance metric. Let f G {Dfi) and let ||/||_h-^» < p with 
0 < a < 1. IfT>f 4 > is zero on for every G [0,27r], then there exists a constant c (a) 
such that: 

||/||l 2 (r 2 ) < c(a) (36) 
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Proof. Let / be defined as in Theorem and let | ■ be the seminorm defined in 
Lemma Let a + 3/2 = m + a for some integer m and 0 < a < 1. Then we have: 

II’^/IIl2([-1,1]xS1) = [ ll’^/0llL2([-l,l])d0 

Js^ 

< / ff 



51 J J[_l,l]x[-l,l] 



dxd?/d0 


(37) 


<ci {afh^-mWl 
< C2 [af 

for 0 < a < 1 with ||/||_fft« < p. Applying Theorem]^ we have: 

3/2 _ 

||/||l2(k 2) < C3 (a)p“+'/"P/lll27[-i,i]xsi) 

< c (a) h°'p 

which completes the proof. 


(38) 


□ 


This last result tells us that given a finite set of detectors with a disc diameter 


sampling determined by equation (31) and with h being a measure of the uniformity 


of the sample, the least squares error in our solution is bounded above by c (a) h°‘p 
with the a-priori knowledge that ||/||j 7 c, < p for some 0 < a < 1. 


3 The physical model 

In this section we present an accurate physical model and a possible approximate 
model which allows for the proposed reconstruction method. We consider an intensity 
of photons scattering from a point u as illustrated in figure The intensity of photons 
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V 



Figure 5: A scattering event with initial photon energy E\ from a source s scattered 
to d with energy Eg. The dashed line displays the original path of the photon to a 
detector v. 


scattered from u to d with energy Eg is: 

I (u, d, Eg) = Iq {Ex) exp {u) dV 

da 


X ^ (Eg^uj) S (q) exp ( - / ) dfi„,d 


(39) 


where Jq is the initial intensity, which depends on the energy Ex (see hgure for 
an example polychromatic spectrum), is the linear attenuation coefficient, which 
is dependant on the energy E and the atomic number of the target material. Here 
Ue {u) dV is the number of electrons in a volume dV around the scattering point u. So 
He (number of electrons per unit volume) is the quantity to be reconstructed. li and 
I 2 are the line segments connecting s to u and u to d respectively. 

The Klein-Nishina differential cross section dcr/dfl, is dehned by: 


da rUEgV (Eg Ex 2 

— (F/ 5 ,a;) =-1 + cos a; 


dn 


2 \ E 


E\ Eg 


(40) 


where tq is the classical electron radius. This predicts the scattering distribution 
for a photon off a free electron at rest. Given that the atomic electrons typically are 
neither free nor at rest, a correction factor is included, namely the incoherent scattering 
function S {q). Here q = ^ sin { 00 /2) is the momentum transferred by a photon with 
initial energy: 

= l-{Eg/Eo){l-cosu) 

scattering at an angle u, where h is Planck’s constant and c is the speed of light. The 
scattering function S also depends on the atomic number Z, so we set Z = Z^vg to 
some average atomic number as an approximation. 
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Figure 6: A typical Tungsten target spectrum with a ISOkV accelerating voltage. 


For Zavg = 45 (Rhodium) we have the expression; 


(1 + 0.458g) 


2.509 


(42) 


To acquire equation (42) we have extended the least squares fit given in [7j to the 
values of S (g) given in jH]. 

The solid angle subtended by u and d is dehned: 

A r ■ n 


= (43) 

dvr |r|'^ 

where r = d — u, A is the detector area and n is the unit vector normal to the detector 
surface. 

Given our machine geometry and proposed reconstruction method, it is difficult 
to include the more accurate model stated above as an additional weighting to our 
integral equations (as in done in [3] for example) while allowing for the same inversion 


formulae. So we average equation (39) over the scattering region = Di ^ADr, for 

T) ’ ' 


each p and <p. Here and Dr are as defined in section where r is fixed depending 
on the machine specifications. 

Let / (u, d, Eg) = I (m; p, 0) = P (n; p, 4>) Ue (u) dV. Here P = P {u; p, 0) depends 
on the scattering point u and p and (j) as defined in section When Rp 7^ 0, P and 
0 determine the detector position d and the measured energy Eg. 

We have: 


-Pavg 


A (Rp^^) J 


P {u;p, 0) du 


(44) 


v-<l> 
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which gives the average of P over Rp,^. Here A (Rp,,ji) denotes the area of Rp,</,. Let 
/ : —)■ M be an example density with snpport contained in Dj., and let: 


h = sC 


P {u;p, 0) du 


(45) 


R. 


■p.4, 


be the scattered intensity measnred for a constant density C over Rp,,/,, where s is 
the (constant) slice thickness. Then if we assnme that the scattering probability is 
constant and eqnal to Pavg(p, 4>) thronghont each scattering region Rp.,/,, the absolnte 
error in onr approximation would satisfy: 


\Im{p, 0) - sPavg {P, (/>) Vlf{p, 0)1 < 


h ( max /(a;),p,0 - h min /(a;),p,0 


6 ) 


for all (p, 0) G Z~^. Here : Z~^ —)■ M is the intensity of photons we measure. So 
provided that the range of the density values is small over the majority of scattering 
regions considered, the averaged model given above will have a similar level of accuracy 


to the more precise model given in equation (39). 


If the linear attenuation coefficient p is known a-priori, then the exponential terms 


of equation (39) may be included in P. Otherwise we may approximate: 


exp ( - y exp j pe 


exp - / Pe, 


(47) 


where ly is the line segment from s to the detector in the forward direction v (see 
hgurej^. This is the approximation made in m- By the Beer-Lambert law, we have: 

-t f \ (48) 


/.(£.) “'"PI-"*"' 

where ly {Ex) is the recorded straight through intensity. 

To account for the physical modelling, we would divide the data by sPavg to cal¬ 
culate approximate values for Pi/ and hence for P/. 


3.1 Determining the atomic number 

With the proposed machine configuration, we can show that the data collected in 
the light field determines the linear attenuation coefficient p uniquely (this is the 
standard 2D reconstruction problem). With the additional information provided by 
our theory, we show under the right assumptions that the atomic number of the target 
is determined uniquely by the full data (light plus dark held). 
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The electron density rie and the linear attennation coefficient /r are related via the 
formnla: 

{E, Z) = Ueae {E, Z) (49) 

where cxe is the total cross section per electron. The cross section cxe is continnons and 
monotone increasing as a fnnction of Z on [1, Zmax], where Z m»^ oc -^/E. For example, 
when E = lOOkeV we can calcnlate .^max = 86. With this in mind, we £t a smooth 
cnrve to known values of cXe given in the atomic data tables [9]. This allows us to 
calculate values of for non integer Z. See figure 



Figure 7: We have presented our fit for erg for E = lOOkeV up to Z = 86 where a 
sudden dip in the cXe values occurs. The tabulated values of ae given in j9| are shown 
alongside the fitted curve. 

The formula for the fit presented for a fixed energy E = lOOkeV is: 

ae (Z) = ap (Z)+a, (Z) = (l.51 X 10“®Z^'^2-0.221ogZ^^('g_4g ^ ^ ^ ^g-4 _ ^-0.50j ^1.57j 

(50) 

This was obtained via a combination of the formula for a^ (the total scattering cross 
section) presented by Jackson and Hawkes in (TO] and the suggested fit for a^ (the 
photoelectric cross section) given in [9j. Fits for energies other than E = lOOkeV are 
also possible via the same fitting method. 
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From our theory we know that the data determines he and He uniquely, where 
E < -Emax- If we assume that atomic numbers Z > are not present in the target 
material, then it is clear from the above arguments that the atomic number of the 
target is uniquely determined. Without this assumption the atomic number would be 
limited to a range of values. If we reconstruct both ^e and Ug for a suitably high 
energy E, we can then calculate values for Z from our curve fit for We will test 
this additional method in our results also. 


4 Results 


To test our reconstruction methods, let us consider the water bottle cross section / 
and the corresponding function / represented in hgure We calculate values of Vf 
for p in the range [0,1] and for 0 G {jfo’''' ’ These were calculated using the 
exact formula for the area of intersection of two discs. We approximate the derivative 
of Vf with respect to p as the finite difference: 


^ A.\ T>/ (p + h, 0) - T>/ (p, 0) 

ifp^f =- h - 


(51) 


for a chosen step size h. To reconstruct / we apply the Matlab function “iradon”, which 
hlters (choosing from a selection of filters pre-coded by Matlab) and backprojects the 
projection data Rf = —-^Vf to recover /. We then make the necessary change in 
coordinates to produce our density image /. In the absence of noise we hnd our results 
to be satisfactory. See hgure 

Let us now perturb the calculated values of Vf slightly such as to simulate random 
noise. We multiply each exact value of Rf by a pseudo random number in the range 


[1 — + q§y] (we use the C-|--|- function “rand” to generate random numbers), 

where %err is the desired amount of percentage error to be added. In this case, even 
with a relatively small amount of added noise, the data must be smoothed sufficiently 


before applying approximation (51). To smooth the data, we apply a simple moving 
average hlter and calculate any intermediate values via a shape preserving cubic in¬ 
terpolation method (“pchip” interpolation in Matlab). We expect this interpolation 
method to preserve the monotonicity of the data (monotone decreasing) as a function 
of p. To illustrate this technique we refer to hgure [OT We have presented our recon- 
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structions after smoothing with 2%, 10% and 50% added noise in hgures 11, 12 and 


Here, we have reconstructed / from a single view point, using data collected from 
a single source projection. With the proposed machine configuration however, there 
are a number of views from which / may be reconstructed. So we take an average 
over 360 views (for source positions at equal vr/lSO intervals over the range [0,27r]). 


Our results are presented in hgures and Here we see an improvement in the 
signal-to-noise-ratio. The rotational symmetry of / about the centre of the circular 
region of /’s support is also recovered. 

Let /avg be the average of the non zero pixel values shown in the left hand image 


of hgure M and let Z = 7.420 be the effective atomic number for water. Then we can 
calculate /avg ~ 1.033 and using equation (49) we can calculate the total cross section 
to be: 

„iF. 7 49n'l 

(52) 


cTe {E, 7.420) = = 0.493 


1.033 

for E = lOOkeV assuming no additional error. Based on our curve ht for cTe(100, Z\ 
this would yield a reconstructed atomic number of Z = 0.886, which differs from the 
accepted value by 88%. For the remaining averaged density reconstructions the /avg 
and Z values are given in the hgure caption. 

We have presented reconstructions of a density which is homogeneous where it 
is not known to be zero. To give an inhomogeneous example, we have presented 
reconstructions with varying levels of added noise of a simulated hollow tube cross 


section in hgures [T^ and 17 


We can summarize our method as follows: 

1. Measure the scattered intensity energy E^ and divide by Pavg and the slice thick¬ 
ness to calculate values for Vf. 


2. Smooth the data sufficiently and apply approximation (51) to calculate values 
for Rf. 

3. Reconstruct / by hltered backprojection and recover / from the dehnition given 
in Theorem [H 

4. Average over a number of source views to improve the image quality and set / 
to 0 outside its support. 
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Figure 8: A water bottle cross section / is simulated as a circular region of uniform 
density 1 on the left. The function / as defined in Theorem is shown on the right. 



Figure 9: A reconstruction of / in the absence of added noise is shown on the left. 
We have backprojected from 180 views with the default Ram-Lak cropped filter. The 
corresponding pixel values of / are presented on the right. Both / and / are set to 0 
outside of their support. 
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Figure 10: On the left we have plotted values of Vf [p, 0) for p > 0 with 10% random 
noise added. On the right we have applied a simple moving average filter to the 
simulated data and taken a subsample of the smoothed data before interpolating as 
specihed earlier. The exact values are presented alongside the fitted values in the right 
hand figure. 



50 100 150 200 100 200 300 400 500 BOO 


Figure 11: On the left we have a reconstruction of / after smoothing with 2% added 
noise. We have again backprojected from 180 views, although here we have multiplied 
the standard ramp filter by a Hamming window to reduce the high frequency noise. 
The corresponding pixel values for / are presented on the right. 
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Figure 12: On the left, a reconstruction of / after smoothing with 10% added noise. 
We have multiplied the ramp hlter by a Hamming window and backprojected from 
180 views. The corresponding pixel values for / are displayed on the right. 



Figure 13: A reconstruction of / after smoothing with 50% added noise is shown on 
the left. We have multiplied the ramp hlter by a Hamming window and backprojected 
from 180 views. The corresponding pixel values for / are displayed on the right. 



























































5 Conclusion 


We have proposed a new fast method to determine the electron density in x-ray scan¬ 
ning applications, with a hxed energy sensitive detector machine conhguration where 
it is possible to measnre photon intensity in the dark held. We have shown that 
the density may be reconstructed analytically using the Compton scattered intensity. 
This method does not require the photon source to be monochromatic as is the case 
in recent literature, which is important from a practical standpoint as it may not be 
reasonable to assume a monochromatic source in some applications. Also if the source 
is monochromatic we cannot gain any insight into the energy dependence of the at¬ 
tenuation coefficient, which would rule out the recent advances in image rendering 
HE], where a combination of multivariate and cluster analysis can be used to render 
a colour x-ray image. 

Using Sobolev space estimates, we have determined an upper bound for the least 
squares error in our solution in terms of the least squares error in our data. This work 
is based on the approach taken by Natterer in [T^ . 

We have shown, under the right assumptions, that the atomic number of the target 
is determined uniquely by the full data. With this theory in place we intend to pursue 
a more practical means to reconstruct the atomic number Z, as the graph reading 
method used in the present paper was ineffective in giving an accurate reconstruction 
for Z. 

We summarize our method to recover the density image in section and we recon¬ 
struct a simulated water bottle cross section via a possible practical implementation 
of this method. In this simple case the smoothing method (simple moving average) 
applied was effective and we were able to reconstruct a circular cross section of ap¬ 
proximately uniform density. Although in the presence of noise the pixel values of our 
reconstructed density image on average differed from the original values by as much as 
15%. We have also provided reconstructions of a simulated hollow tube cross section. 
In this case the inner edge of the tube cross section appeared quite blurred in the 
reconstruction when noise was added to the simulated data. We performed a num¬ 
ber of trial reconstructions with different randomly generated datasets. The results 
presented in this paper are typical of our trial results. 
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We hope also to test our methods through experiment. For example, if we were 
to take an existing x-ray machine of a similar conhguration to that discussed in the 
present paper, and attach energy sensitive detectors alongside the existing detectors or 
if we were to replace them, then we could see how closely our forward problem models 
the intensity of photons measured in the dark field in practice. 
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Appendix — The RTT80; An example application in 
threat detection 


The RTT80 (real time tomography) X-ray scanner is a switched source, offset detector 
CT machine designed with the aim to scan objects in real time. Developed by Rapiscan 
systems, the RTT80 is currently used in airport security screening of baggage. 

The RTT80 consists of a single fixed ring of polychromatic X-ray sources and 
multiple offset rings of detectors, with a conveyor belt and scanning tunnel (within 
which the scanned object would be placed) passing through the centre of both sets of 
rings. See figure [T^ If the detectors are energy sensitive, then in this case we have the 
problem of reconstructing a density slice supported within the scanning tunnel from 
its integrals over toric sections, with tips at the source and detector locations. We 
wish to check whether it is reasonable to approximate a set of toric section integrals 
as integrals over discs whose boundaries intersect a given source point, as then we can 
apply our proposed reconstruction method to reconstruct the density slice analytically. 

Let us refer to figure and let Dp be defined as in section We define the toric 
sections = D^^nDp,^, Tl^ = D^,^UDp,^ and = D^^UDp,^. 

Let A(5') denote the area of a set S' C and let T C denote the set of points 
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within our ROI (region of interest, i.e the scanning tunnel). For a large sample of 
discs, we will check for every disc in the sample, whether 3i G {1,2, 3,4} such 
that A{Dp^^ n T) ^ ^ 

Let Dr be dehned as in section Then if we consider the machine specihcations 
for the RTT80, we can calculate r = 6.75 and the difference in radius between the 
detector ring and the scanning tunnel to be 0.375. See hgure For our test, we 
consider a sample of 36000 discs with diameters p = 1.375 + for 1 < z < 100 and 
0 = for 1 < j < 360. We have chosen p G [1.375, 6.375] and 0 G [0, 27r] values in 
a range sufficient to determine a unique density slice image for densities supported on 
T. Refer to Corollary [Tj For each of our chosen p and 0 value pairs, the difference: 


mm (A(D,,t n T) - n T)) « 10-‘« (53) 

was found to be negligible. Let / : —)■ M be an example density slice with support 

contained in T. Then for any disc Dp^^ in our sample we have: 


/ = 


/ = 


/ = 


/ 


(54) 


'D, 


v,<l> 




IT’- ,r\T 

p,<p 


pA 


which holds for some i G (1,2, 3,4}. So, the integral of / over ^ is equal to at 
least one of four toric section integrals over /. Assuming also that there is little error 
implied by our physical approximations (these are discussed in detail in section]^. 


the integral (54) would be determined approximately by at least one of four data sets, 
namely the photon intensity measured for two possible energy levels at two possible 
detector locations (d^ ^ or dp ,p). Thus, given that the inverse disc transform is only 


mildly ill posed (this was determined to be the case in section 2.1, based on the criteria 
given by Natterer in [I2]), it seems that we should be able obtain a satisfactory density 
image reconstruction in this application. 

In airport baggage screening, we are interested in identifying a given material as 
either a threat or non-threat. Let Ue be the electron density and let Z denote the 
effective atomic number. We dehne the threat space to be the set of materials with 
(ue, Z) G T, where T C [0, cx)) x [1,100] is the class of threat (ue, Z) pairs. For a given 
suspect material, we can apply the methods presented in this paper to reconstruct 
Ue and Z. Then if (ug, Z) G T, we can identify the suspect material as a potential 
threat. We note that although we failed to obtain an accurate Z reconstruction in the 
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present paper, we aim to show that a more precise determination of Z is possible in 
future work. Also, the reconstruction methods we have presented should be fast to 
implement as they are largely based on the hltered back-projection algorithm. This 
is important in an application such as airport baggage screening, as we require the 
threat detection method we apply to not only be accurate in threat identihcation, but 
to also be an efficient process. 
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Figure 18: The RTT80 machine conhguration is displayed. The source-detector ring 
offset is relatively small and will be modelled as zero. The RTT80’s relative dimensions 
(the source ring, detector ring and scanning tunnel radii) are presented to the left of 
the diagram. 



Figure 19: The RTT80 configuration is displayed. The origin is denoted by O as in 
section This is where a source is located. A disc with boundary intersecting O 
and two detector points ^ and is shown to have a non empty intersection with 
the set of points in our ROI (the scanning tunnel). The disc is the reflection of 
Dp^^ in the line segment connecting O to Similarly for 
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